#########################################################################################################
# Replication file for Lehrer, Juhl, Stoeckle:                                                          #
# "Assessing the relative influence of party unity on vote choice: Evidence from a conjoint experiment" #
#                                                                                                       #
# Created by Roni Lehrer                                                                                #
# Last edited on 1 August 2022                                                                          #
#########################################################################################################

# cleaning directory
rm(list=ls())

#loading required packages
library(data.table)
library(survival)
library(MASS)
library(tidyverse)
library(knitr)
library(kableExtra)
library(cjoint)

#if tinytex is not installed yet, please install run:

install.packages(tinytex)
tinytex::install_tinytex()
library(tinytex)

#load data
load("Leher_replication_data.Rdata")

# Setting up log
log_file <- file("Leher_replication_log.txt") # Name of log file
sink(log_file, append = TRUE, type = "output") # Writing console output to log file
sink(log_file, append = TRUE, type = "message") # Writing messages to log file


#set a seed
set.seed(363048) #from random.org

#############################
# Baseline Regression Model #
#############################
model <- clogit(chosen~
                  dist +
                  critique +
                  parliament +
                  conference +
                  reform +
                  role +
                  gender +
                  age +
                  job +
                  strata(id_screen) +
                  cluster(id_g),
                data=dat, method="efron", robust=TRUE)

# draw from sampling distribution (required for several Figures below)
betas <- mvrnorm(1000, coef(model), vcov(model)) #the vcov is already robust

######################
# Figure 1: AME Plot #
######################
compute.pdata <- function(model) {
    #get data for AMCEs plot in tidy format
    coefs <- coef(model) #get cofficients
    ses <- summary(model)$coefficients[,4] #get ROBUST SEs
    names <- rownames(summary(model)$coefficients) #get names of coefs
    pdata <- data.table(mean=(1/(1+exp(-coefs)))-.5, #compute probability difference to all 0's scenario
                        lower=(1/(1+exp(-(coefs-(1.96*ses)))))-.5, #compute corresponding CI
                        upper=(1/(1+exp(-(coefs+(1.96*ses)))))-.5, #compute corresponding CI
                        names=names
                        )
    return(pdata)
}

pdata <- compute.pdata(model = model)

#set nicer names (as they appear in Figure 1)
pdata$nice.names <- ifelse(pdata$names=="dist1"                               , "1", "")
pdata$nice.names <- ifelse(pdata$names=="dist2"                               , "2", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="dist3"                               , "3", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="dist4"                               , "4", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="critiqueFormer party leader"   , "Former party leader", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="critiqueParty factions"               , "Party factions", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="critiqueRank-and-file members", "Rank-and-file members", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="parliamentDivided voting"            , "Divided voting", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="conferenceNeither united nor divided", "Neither united nor divided", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="conferenceDivided"                   , "Divided", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="reformLow"                           , "Low", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="roleOpposition party"            , "Opposition party", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="roleJunior coalition partner", "Junior coalition partner",
                           pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="genderMale"                          , "Male", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="age56 years"                         , "56 years", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="age74 years"                         , "74 years", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="jobActivist"                         , "Activist", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="jobLawyer"                           , "Lawyer", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="jobPolitician"                       , "Politician", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="jobEntrepreneur"                     , "Entrepreneur", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="jobEmployee (retired)"          , "Employee (retired)", pdata$nice.names)

#addd new lines for reference categories and headers
pdata <- rbindlist(list(pdata,
                   data.table(mean=0, lower=0, upper=0, names="dist0", nice.names="0"),
                   data.table(mean=0, lower=0, upper=0, names="critiqueNone", nice.names="None"),
                   data.table(mean=0, lower=0, upper=0, names="parliamentUnited", nice.names="United voting"),
                   data.table(mean=0, lower=0, upper=0, names="conferenceUnited", nice.names="United"),
                   data.table(mean=0, lower=0, upper=0, names="reformHigh", nice.names="High"),
                   data.table(mean=0, lower=0, upper=0, names="rolePM party", nice.names="PM party"),
                   data.table(mean=0, lower=0, upper=0, names="genderFemale", nice.names="Female"),
                   data.table(mean=0, lower=0, upper=0, names="jobEmployee", nice.names="Employee"),
                   data.table(mean=0, lower=0, upper=0, names="age38 years", nice.names="38 years"),
                   data.table(mean=0, lower=0, upper=0, names="header", nice.names="Ideological Distance"),
                   data.table(mean=0, lower=0, upper=0, names="header", nice.names="Party Role"),
                   data.table(mean=0, lower=0, upper=0, names="header", nice.names="Intra-Party Critique"),
                  data.table(mean=0, lower=0, upper=0, names="header", nice.names="Parliamentary Voting Behavior"),
                   data.table(mean=0, lower=0, upper=0, names="header", nice.names="Behavior at Party Congress"),
                   data.table(mean=0, lower=0, upper=0, names="header", nice.names="Reform Clarity"),
                   data.table(mean=0, lower=0, upper=0, names="header", nice.names="Candidate: Gender"),
                   data.table(mean=0, lower=0, upper=0, names="header", nice.names="Candidate: Age"),
                   data.table(mean=0, lower=0, upper=0, names="header", nice.names="Candidate: Occupation")))

## # check new names versus those from the analysis
## pdata[,list(names, nice.names)]

#assign correct y values to each line
pdata$y <- ifelse(pdata$names=="jobActivist", 1, NA)
pdata$y <- ifelse(pdata$names=="jobEmployee (retired)", 2, pdata$y)
pdata$y <- ifelse(pdata$names=="jobPolitician", 3, pdata$y)
pdata$y <- ifelse(pdata$names=="jobLawyer", 4, pdata$y)
pdata$y <- ifelse(pdata$names=="jobEntrepreneur", 5, pdata$y)
pdata$y <- ifelse(pdata$names=="jobEmployee", 6, pdata$y)
pdata$y <- ifelse(pdata$nice.names=="Candidate: Occupation", 7, pdata$y)
pdata$y <- ifelse(pdata$names=="age74 years", 8, pdata$y)
pdata$y <- ifelse(pdata$names=="age56 years", 9, pdata$y)
pdata$y <- ifelse(pdata$names=="age38 years", 10, pdata$y)
pdata$y <- ifelse(pdata$nice.names=="Candidate: Age", 11, pdata$y)
pdata$y <- ifelse(pdata$names=="genderMale", 12, pdata$y)
pdata$y <- ifelse(pdata$names=="genderFemale", 13, pdata$y)
pdata$y <- ifelse(pdata$nice.names=="Candidate: Gender", 14, pdata$y)
pdata$y <- ifelse(pdata$names=="roleOpposition party", 15, pdata$y)
pdata$y <- ifelse(pdata$names=="roleJunior coalition partner", 16, pdata$y)
pdata$y <- ifelse(pdata$names=="rolePM party", 17, pdata$y)
pdata$y <- ifelse(pdata$nice.names=="Party Role", 18, pdata$y)
pdata$y <- ifelse(pdata$names=="reformLow", 19, pdata$y)
pdata$y <- ifelse(pdata$names=="reformHigh", 20, pdata$y)
pdata$y <- ifelse(pdata$nice.names=="Reform Clarity", 21, pdata$y)
pdata$y <- ifelse(pdata$names=="conferenceDivided", 22, pdata$y)
pdata$y <- ifelse(pdata$names=="conferenceNeither united nor divided", 23, pdata$y)
pdata$y <- ifelse(pdata$names=="conferenceUnited", 24, pdata$y)
pdata$y <- ifelse(pdata$nice.names=="Behavior at Party Congress", 25, pdata$y)
pdata$y <- ifelse(pdata$names=="parliamentDivided voting", 26, pdata$y)
pdata$y <- ifelse(pdata$names=="parliamentUnited", 27, pdata$y)
pdata$y <- ifelse(pdata$nice.names=="Parliamentary Voting Behavior", 28, pdata$y)
pdata$y <- ifelse(pdata$names=="critiqueParty factions", 29, pdata$y)
pdata$y <- ifelse(pdata$names=="critiqueFormer party leader", 30, pdata$y)
pdata$y <- ifelse(pdata$names=="critiqueRank-and-file members", 31, pdata$y)
pdata$y <- ifelse(pdata$names=="critiqueNone", 32, pdata$y)
pdata$y <- ifelse(pdata$nice.names=="Intra-Party Critique", 33, pdata$y)
pdata$y <- ifelse(pdata$names=="dist4", 34, pdata$y)
pdata$y <- ifelse(pdata$names=="dist3", 35, pdata$y)
pdata$y <- ifelse(pdata$names=="dist2", 36, pdata$y)
pdata$y <- ifelse(pdata$names=="dist1", 37, pdata$y)
pdata$y <- ifelse(pdata$names=="dist0", 38, pdata$y)
pdata$y <- ifelse(pdata$nice.names=="Ideological Distance", 39, pdata$y)

#omit headers from plotting coefficients
pdata$mean <- ifelse(pdata$names=="header", NA, pdata$mean)
pdata$lower <- ifelse(pdata$names=="header", NA, pdata$lower)
pdata$upper <- ifelse(pdata$names=="header", NA, pdata$upper)

p <- ggplot(data=pdata, group=y)
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=.5, ymax=6.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=7.5, ymax=10.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=11.5, ymax=13.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=14.5, ymax=17.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=18.5, ymax=20.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=21.5, ymax=24.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=25.5, ymax=27.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=28.5, ymax=32.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=33.5, ymax=38.5), fill="grey")
p <- p + geom_segment(aes(x=0, xend=0, y=0, yend=38.5))
p <- p + geom_vline(xintercept=-.1, color="white")
p <- p + geom_vline(xintercept=-.2, color="white")
p <- p + geom_vline(xintercept=-.3, color="white")
p <- p + geom_vline(xintercept=-.4, color="white")
p <- p + geom_vline(xintercept=-.5, color="white")
p <- p + geom_segment(aes(x=lower, xend=upper, y=y, yend=y))
p <- p + geom_point(aes(x=mean, y=y))
p <- p + scale_y_continuous("", position="right",
                            breaks=pdata[names!="header", y],
                            labels=pdata[names!="header", nice.names])
p <- p + theme_bw()
p <- p + theme(panel.grid.major = element_blank(),
               panel.grid.minor = element_blank(),
               panel.grid.major.x = element_line(color="white"),
               panel.grid.minor.x = element_line(color="white"),
               axis.line.x = element_line(),
               panel.border = element_blank(),
               text = element_text(size=14))
p <- p + scale_x_continuous("Average Marginal Component Effect",
                            breaks=seq(0, -.5, -.1))
p <- p + annotate("text", y=pdata[names=="header", y],
                                   x=-.25,
                  label=pdata[names=="header", nice.names])
p <- p + coord_cartesian(xlim = c(-.49, 0),
                         ylim = c(1.5, 37.6),
                         clip = "on")
p
ggsave("Figure1.pdf", units="cm", width=17.8, height=17.8)

##########################
# Figure 2 and Figure A5 #
##########################
compute_probs_one_vs_second <- function(x) {
#Condition Logit: Compute probability that option A or B are chosen
  prob.first <- exp(x[1])/(sum(exp(x)))
  return(c(prob.first, 1-prob.first))
}

#compute predicted probabilities for all scenarios (Fig 2 and Fig A5)
#set scenarios
scenarios <- data.table(dist=c("0", "1", "1", "1", "2"),
              critique= c("Party factions", "Party factions", "Party factions", "None", "None"),
              parliament=c("Divided voting", "Divided voting", "Divided voting", "United voting",
                           "United voting"),
              conference= c("Divided",  "Divided", "United", "Neither united nor divided", "United"),
                      reform=c("High"),
                      role=c("PM party"),
                      gender=c("Female"),
                      age=c("38 years"),
                      job=c( "Employee"),
                      id_screen=1)

# create design matrix that has indicator vars instead of factors
scenarios.expanded <- model.matrix(model, data=scenarios)

# compute linear predictor
linpred <- scenarios.expanded%*%t(betas)

# find median probability and CIs for plotting
scenarios$prob <- apply(linpred, 1, median)
scenarios$upper <- apply(linpred, 1, quantile, p=.975)
scenarios$lower <- apply(linpred, 1, quantile, p=.025)

scenarios$party <- c("1", "2", "3", "4", "5") #assign names from Table 2

#transform linear predictor into probabilities (including CIs) for each pair
predprob1_2 <- apply(scenarios[party%in%c("1", "2"), .(prob, lower, upper)], 2, compute_probs_one_vs_second)
predprob1_3 <- apply(scenarios[party%in%c("1", "3"), .(prob, lower, upper)], 2, compute_probs_one_vs_second)
predprob1_4 <- apply(scenarios[party%in%c("1", "4"), .(prob, lower, upper)], 2, compute_probs_one_vs_second)
predprob1_5 <- apply(scenarios[party%in%c("1", "5"), .(prob, lower, upper)], 2, compute_probs_one_vs_second)

#The following is for Fig 2, for Fig A5 (Party 1 vs. Party 5) see below.
pdata <- data.table(prob=c(predprob1_2[,1], predprob1_3[,1], predprob1_4[,1]),
                    lower=c(predprob1_2[,2], predprob1_3[,2], predprob1_4[,2]),
                    upper=c(predprob1_2[,3], predprob1_3[,3], predprob1_4[,3]),
                    party=c("1", "2", "1", "3", "1", "4"),
                    case=c("1vs2", "1vs2", "1vs3", "1vs3", "1vs4", "1vs4")
                    )
pdata$y <- ifelse(pdata$case=="1vs2", 3,
           ifelse(pdata$case=="1vs3", 2, 1))

p <- ggplot(data=pdata[party=="1"], aes(group=case))
p <- p + theme_bw()
p <- p + geom_rect(aes(xmin=-1, xmax=101, ymin=.75, ymax=1.23), fill="grey")
p <- p + geom_rect(aes(xmin=-1, xmax=101, ymin=1.75, ymax=2.23), fill="grey")
p <- p + geom_rect(aes(xmin=-1, xmax=101, ymin=2.75, ymax=3.23), fill="grey")
p <- p + theme(panel.grid.major = element_blank(),
               panel.grid.minor = element_blank(),
               panel.border = element_blank(),
               axis.text = element_blank(),
               axis.ticks = element_blank())
p <- p + geom_segment(aes(x=0, xend=0, y=.75, yend=1.25), linetype="dotted")
p <- p + geom_segment(aes(x=0, xend=0, y=1.75, yend=2.25), linetype="dotted")
p <- p + geom_segment(aes(x=0, xend=0, y=2.75, yend=3.25), linetype="dotted")
p <- p + geom_segment(aes(x=25, xend=25, y=.75, yend=1.25), linetype="dotted")
p <- p + geom_segment(aes(x=25, xend=25, y=1.75, yend=2.25), linetype="dotted")
p <- p + geom_segment(aes(x=25, xend=25, y=2.75, yend=3.25), linetype="dotted")
p <- p + geom_segment(aes(x=50, xend=50, y=.75, yend=1.25), linetype="dotted")
p <- p + geom_segment(aes(x=50, xend=50, y=1.75, yend=2.25), linetype="dotted")
p <- p + geom_segment(aes(x=50, xend=50, y=2.75, yend=3.25), linetype="dotted")
p <- p + geom_segment(aes(x=75, xend=75, y=.75, yend=1.25), linetype="dotted")
p <- p + geom_segment(aes(x=75, xend=75, y=1.75, yend=2.25), linetype="dotted")
p <- p + geom_segment(aes(x=75, xend=75, y=2.75, yend=3.25), linetype="dotted")
p <- p + geom_segment(aes(x=100, xend=100, y=.75, yend=1.25), linetype="dotted")
p <- p + geom_segment(aes(x=100, xend=100, y=1.75, yend=2.25), linetype="dotted")
p <- p + geom_segment(aes(x=100, xend=100, y=2.75, yend=3.25), linetype="dotted")
p <- p + annotate("text", x=0, y=.7, label="100:0")
p <- p + annotate("text", x=0, y=1.7, label="100:0")
p <- p + annotate("text", x=0, y=2.7, label="100:0")
p <- p + annotate("text", x=25, y=.7, label="75:25")
p <- p + annotate("text", x=25, y=1.7, label="75:25")
p <- p + annotate("text", x=25, y=2.7, label="75:25")
p <- p + annotate("text", x=50, y=.7, label="50:50")
p <- p + annotate("text", x=50, y=1.7, label="50:50")
p <- p + annotate("text", x=50, y=2.7, label="50:50")
p <- p + annotate("text", x=75, y=.7, label="25:75")
p <- p + annotate("text", x=75, y=1.7, label="25:75")
p <- p + annotate("text", x=75, y=2.7, label="25:75")
p <- p + annotate("text", x=100, y=.7, label="0:100")
p <- p + annotate("text", x=100, y=1.7, label="0:100")
p <- p + annotate("text", x=100, y=2.7, label="0:100")
p <- p + annotate("text", x=5, y=3.4, label="Party 1", fontface = "bold")
p <- p + annotate("text", x=95, y=3.4, label="Party 2", fontface = "bold")
p <- p + annotate("text", x=5, y=2.4, label="Party 1", fontface = "bold")
p <- p + annotate("text", x=95, y=2.4, label="Party 3", fontface = "bold")
p <- p + annotate("text", x=5, y=1.4, label="Party 1", fontface = "bold")
p <- p + annotate("text", x=95, y=1.4, label="Party 4", fontface = "bold")
p <- p + ylab("")
p <- p + xlab("")
p <- p + geom_rect(aes(xmin=100-(lower*100), xmax=100-(upper*100), ymin=y-.3, ymax=y+.3), fill="grey37")
p <- p + geom_segment(aes(x=100-(prob*100), xend=100-(prob*100), y=y-.32, yend=y+.32))
p <- p + annotate("text", y=3.4, x=39.7, label="60:40")
p <- p + annotate("text", y=2.4, x=55, label="45:55")
p <- p + annotate("text", y=1.4, x=58, label="42:58")
p
ggsave(file="Figure2.pdf", units="cm", height=10, width=17)

#The following is for Fig A5:
pdata <- data.table(prob=predprob1_5[,1],
                    lower=predprob1_5[,2],
                    upper=predprob1_5[,3],
                    party=c("1", "5"),
                    case=c("1vs5", "1vs5")
                    )
pdata$y <- 1

p <- ggplot(data=pdata[party=="1"], aes(group=case))
p <- p + theme_bw()
p <- p + geom_rect(aes(xmin=-1, xmax=101, ymin=.75, ymax=1.23), fill="grey")
p <- p + theme(panel.grid.major = element_blank(),
               panel.grid.minor = element_blank(),
               panel.border = element_blank(),
               axis.text = element_blank(),
               axis.ticks = element_blank())
p <- p + geom_segment(aes(x=0, xend=0, y=.75, yend=1.25), linetype="dotted")
p <- p + geom_segment(aes(x=25, xend=25, y=.75, yend=1.25), linetype="dotted")
p <- p + geom_segment(aes(x=50, xend=50, y=.75, yend=1.25), linetype="dotted")
p <- p + geom_segment(aes(x=75, xend=75, y=.75, yend=1.25), linetype="dotted")
p <- p + geom_segment(aes(x=100, xend=100, y=.75, yend=1.25), linetype="dotted")
p <- p + annotate("text", x=0, y=.7, label="100:0")
p <- p + annotate("text", x=25, y=.7, label="75:25")
p <- p + annotate("text", x=50, y=.7, label="50:50")
p <- p + annotate("text", x=75, y=.7, label="25:75")
p <- p + annotate("text", x=100, y=.7, label="0:100")
p <- p + annotate("text", x=5, y=1.4, label="Party 1", fontface = "bold")
p <- p + annotate("text", x=95, y=1.4, label="Party 6", fontface = "bold")
p <- p + ylab("")
p <- p + xlab("")
p <- p + geom_rect(aes(xmin=100-(lower*100), xmax=100-(upper*100), ymin=y-.3, ymax=y+.3), fill="grey37")
p <- p + geom_segment(aes(x=100-(prob*100), xend=100-(prob*100), y=y-.32, yend=y+.32))
p <- p + annotate("text", y=1.4, x=37, label="63:37")
p
ggsave(file="FigureA5.pdf", units="cm", height=5, width=17)

#####################################
# Table A3: Regression Coefficients #
#####################################
coefs <- coef(model)
rses <- as.character(round(summary(model)$coefficients[,4], 3))
zs <- coef(model)/summary(model)$coefficients[,4]
ps <- 2*pnorm(-abs(zs))
names <- rownames(summary(model)$coefficients)

reg_table <- data.table(names, coefs, rses, zs, ps)

reg_table$coefs_stars <-  ifelse(reg_table$ps<0.05,
                          ifelse(reg_table$ps<0.01,
                          ifelse(reg_table$ps<0.001, 
                          paste0(sprintf("%.3f", reg_table$coefs), "***"),
                          paste0(sprintf("%.3f", reg_table$coefs), "**")),
                          paste0(sprintf("%.3f", reg_table$coefs), "*")),
                                 sprintf("%.3f", reg_table$coefs))

out <- reg_table[,list(names, coefs_stars, rses)]

# add empty rows to include the baseline categories
out <- rbindlist(list(
    data.table(names = "NA",coefs_stars = "", rses = ""),
    out[1:4],
    data.table(names = "NA",coefs_stars = "", rses = ""),
    out[5:7],
    data.table(names = "NA",coefs_stars = "", rses = ""),
    out[8:8],
    data.table(names = "NA",coefs_stars = "", rses = ""),
    out[9:10],
    data.table(names = "NA",coefs_stars = "", rses = ""),
    out[11:11],
    data.table(names = "NA", coefs_stars = "", rses = ""),
    out[12:13],
    data.table(names = "NA",coefs_stars = "", rses = ""),
    out[14:14],
    data.table(names = "NA", coefs_stars = "", rses = ""),
    out[15:16],
    data.table(names = "NA",coefs_stars = "", rses = ""),
    out[17:21],
    data.table(names = "NA",coefs_stars = "", rses = ""),
    data.table(names = "Log Likelihood",
            coefs_stars = as.character(round(summary(model)$loglik[2], 3)),
            rses = ""),
    data.table(names = "N (observations)",
            coefs_stars = as.character(dat[complete.cases(chosen, dist), .(n_observations = length(id_screen))]),
            rses = ""),
    data.table(names = "N (choices)",
            coefs_stars = as.character(dat[complete.cases(chosen, dist), .(n_tasks = length(unique(id_screen)))]),
            rses = ""),
    data.table(names = "N (respondents)",
            coefs_stars = as.character(dat[complete.cases(chosen, dist), .(n_respondents = length(unique(id_g)))]),
            rses = "")
    ))

# names vector for output
names.out <-  c("Ideological distance (reference: 0)", "1", "2", "3", "4",
         "Intra-Party Critique (reference: None)","Rank-and-file members", "Former party leader", "Party factions",
            "Parliamentary voting (reference: United)", "Divided",
            "Behavior at Party Congress (reference: United)", "Neither united nor divided", "Divided",
            "Reform clarity (reference: High)", "Low",
            "Party role (reference: PM party)", "Opposition party", "Junior coalition partner",
            "Candidate's gender (reference: Female)", "Male",
            "Candidate's age (reference: 38 years)", "56 years", "74 years",
            "Candidate's occupation (reference: Employee)", "Activist", "Lawyer", "Politician", "Entrepreneur", "Employee (retired)",
            "",
            "Log Likelihood",
            "N (observations)", "N (choices)", "N (respondents)")

out$names <- names.out

## # check correct names by row.names
## out
# then remove row-names
row.names(out) <- NULL

tab <- kable(out, format = "latex",
      col.names = c("", "coefficient", "robust s.e."),
      caption = "Estimated coefficients from the conditional logistic model",
      label = "reg-output",
      booktabs = T, linesep = "")
tab <-  add_footnote(tab, "***p<0.001; **p<0.01; *p<0.05", notation = "none")
tab <-  kable_styling(tab, font_size = 11,
                latex_options = c("hold_position"),
                position  = "center")
save_kable(tab, file="TableA3.pdf", keep_tex=TRUE) #save table as tex file

#compile tex file to obtain pdf version
tinytex::xelatex("TableA3.tex")

#########################################################
# Figure A1: comparison of clogit and cjoint regression #
#########################################################
#exclude missings
dat.cbc <- dat[complete.cases(chosen, dist), ]

#set new name for baseline model
model_clogit <- model

# custom function to get data for AMCE plot in consistent format
compute_amce_clogit <- function(model) {
  coefs <- coef(model) #get coefficients
  ses <- summary(model)$coefficients[,4] #get ROBUST ses
  names <- rownames(summary(model)$coefficients) #get names of coefs
  pdata <- data.table(amce=(1/(1+exp(-coefs)))-.5, #compute probability difference to all 0's scenario
                      lower=(1/(1+exp(-(coefs-(1.96*ses)))))-.5, #compute corresponding CI
                      upper=(1/(1+exp(-(coefs+(1.96*ses)))))-.5, #compute corresponding CI
                      names=names,
                      coefs = coefs,
                      ses = ses,
                      specification = "clogit"
  )
  return(pdata)
}


pdata_clogit <- compute_amce_clogit(model = model_clogit)

# approach by Hainmueller et al 2014, implemented in the "cjoint" package
model_cjoint <- amce(chosen~
                    gender
                     +age
                     +job
                     +role
                     +critique
                     +parliament
                     +conference
                     +reform
                     +dist,
                     data = dat.cbc, design = "uniform",
                     respondent.varying = NULL, subset = NULL,
                     respondent.id = "id_g", cluster = TRUE, na.ignore=T,
                     weights = NULL, baselines = NULL)


summary(model_cjoint)

# custom function to get data for AMCE plot in consistent format
compute_amce_cjoint <- function(model) {
  coefs <- summary(model_cjoint)$amce[,3] #get coefficients
  ses <- summary(model_cjoint)$amce[,4] #get ses (not robust here)
  names <- paste0(summary(model_cjoint)$amce$Attribute,summary(model_cjoint)$amce$Level)  #get names of coefs
  pdata <- data.table(amce=coefs, #here, AMCE is just the coefficients
                      lower=coefs-(1.96*ses), #compute corresponding CI
                      upper=coefs+(1.96*ses), #compute corresponding CI
                      names=names,
                      coefs = coefs,
                      ses = ses,
                      specification = "cjoint"
  )
  return(pdata)
}

pdata_cjoint <- compute_amce_cjoint(model = model_cjoint)

write.csv(pdata_cjoint, "pdata_conjoint.csv", row.names=FALSE) #save as CSV

# join in one dataframe
pdata_robustness1 <- rbind(pdata_clogit, pdata_cjoint)

#set nice names
pdata_robustness1$nice.names <- ifelse(pdata_robustness1$names=="dist1" , "1", "")
pdata_robustness1$nice.names <- ifelse(pdata_robustness1$names=="dist2" , "2", pdata_robustness1$nice.names)
pdata_robustness1$nice.names <- ifelse(pdata_robustness1$names=="dist3" , "3", pdata_robustness1$nice.names)
pdata_robustness1$nice.names <- ifelse(pdata_robustness1$names=="dist4" , "4", pdata_robustness1$nice.names)
pdata_robustness1$nice.names <- ifelse(pdata_robustness1$names=="critiqueFormer party leader" ,
                                       "Former party leader", pdata_robustness1$nice.names)
pdata_robustness1$nice.names <- ifelse(pdata_robustness1$names=="critiqueParty factions" ,
                                       "Party factions", pdata_robustness1$nice.names)
pdata_robustness1$nice.names <- ifelse(pdata_robustness1$names=="critiqueRank-and-file members",
                                       "Rank-and-file members", pdata_robustness1$nice.names)
pdata_robustness1$nice.names <- ifelse(pdata_robustness1$names=="parliamentDivided voting",
                                       "Divided voting", pdata_robustness1$nice.names)
pdata_robustness1$nice.names <- ifelse(pdata_robustness1$names=="conferenceNeither united nor divided",
                                       "Neither united nor divided", pdata_robustness1$nice.names)
pdata_robustness1$nice.names <- ifelse(pdata_robustness1$names=="conferenceDivided" ,
                                       "Divided", pdata_robustness1$nice.names)
pdata_robustness1$nice.names <- ifelse(pdata_robustness1$names=="reformLow" , "Low", pdata_robustness1$nice.names)
pdata_robustness1$nice.names <- ifelse(pdata_robustness1$names=="roleOpposition party",
                                       "Opposition party", pdata_robustness1$nice.names)
pdata_robustness1$nice.names <- ifelse(pdata_robustness1$names=="roleJunior coalition partner",
                                       "Junior coalition partner", pdata_robustness1$nice.names)
pdata_robustness1$nice.names <- ifelse(pdata_robustness1$names=="genderMale",
                                       "Male", pdata_robustness1$nice.names)
pdata_robustness1$nice.names <- ifelse(pdata_robustness1$names=="age56 years" ,
                                       "56 years", pdata_robustness1$nice.names)
pdata_robustness1$nice.names <- ifelse(pdata_robustness1$names=="age74 years" ,
                                       "74 years", pdata_robustness1$nice.names)
pdata_robustness1$nice.names <- ifelse(pdata_robustness1$names=="jobActivist" ,
                                       "Activist", pdata_robustness1$nice.names)
pdata_robustness1$nice.names <- ifelse(pdata_robustness1$names=="jobLawyer" ,
                                       "Lawyer", pdata_robustness1$nice.names)
pdata_robustness1$nice.names <- ifelse(pdata_robustness1$names=="jobPolitician" ,
                                       "Politician", pdata_robustness1$nice.names)
pdata_robustness1$nice.names <- ifelse(pdata_robustness1$names=="jobEntrepreneur" ,
                                       "Entrepreneur", pdata_robustness1$nice.names)
pdata_robustness1$nice.names <- ifelse(pdata_robustness1$names=="jobEmployee (retired)",
                                       "Employee (retired)", pdata_robustness1$nice.names)

#addd new lines for reference categories and headers
additional.lines1 <- rbindlist(list(
                   data.table(names="dist0", nice.names="0"),
                   data.table(names="critiqueNone", nice.names="None"),
                   data.table(names="parliamentUnited", nice.names="United voting"),
                   data.table(names="conferenceUnited", nice.names="United"),
                   data.table(names="reformHigh", nice.names="High"),
                   data.table(names="rolePM party", nice.names="PM party"),
                   data.table(names="genderFemale", nice.names="Female"),
                   data.table(names="jobEmployee", nice.names="Employee"),
                   data.table(names="age38 years", nice.names="38 years"),
                   data.table(names="header", nice.names="Ideological Distance"),
                   data.table(names="header", nice.names="Party Role"),
                   data.table(names="header", nice.names="Intra-Party Critique"),
                   data.table(names="header", nice.names="Parliamentary Voting Behavior"),
                   data.table(names="header", nice.names="Behavior at Party Congress"),
                   data.table(names="header", nice.names="Reform Clarity"),
                   data.table(names="header", nice.names="Candidate: Gender"),
                   data.table(names="header", nice.names="Candidate: Age"),
                   data.table(names="header", nice.names="Candidate: Occupation")))
additional.lines1$amce <- ifelse(additional.lines1$names!="header", 0, NA)
additional.lines2 <- additional.lines1
additional.lines1$specification <- "clogit"
additional.lines2$specification <- "cjoint"
pdata_robustness1 <- rbindlist(list(pdata_robustness1, additional.lines1, additional.lines2),
                               use.names=TRUE, fill=TRUE)


## # check new names versus those from the analysis
## pdata_robustness1[,list(names, nice.names)]

#assign correct y values for each line in plot
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="jobActivist", 1, NA)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="jobEmployee (retired)", 2, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="jobPolitician", 3, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="jobLawyer", 4, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="jobEntrepreneur", 5, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="jobEmployee", 6, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$nice.names=="Candidate: Occupation", 7, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="age74 years", 8, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="age56 years", 9, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="age38 years", 10, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$nice.names=="Candidate: Age", 11, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="genderMale", 12, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="genderFemale", 13, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$nice.names=="Candidate: Gender", 14, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="roleOpposition party", 15, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="roleJunior coalition partner", 16, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="rolePM party", 17, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$nice.names=="Party Role", 18, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="reformLow", 19, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="reformHigh", 20, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$nice.names=="Reform Clarity", 21, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="conferenceDivided", 22, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="conferenceNeither united nor divided", 23, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="conferenceUnited", 24, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$nice.names=="Behavior at Party Congress", 25, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="parliamentDivided voting", 26, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="parliamentUnited", 27, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$nice.names=="Parliamentary Voting Behavior", 28, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="critiqueParty factions", 29, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="critiqueFormer party leader", 30, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="critiqueRank-and-file members", 31, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="critiqueNone", 32, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$nice.names=="Intra-Party Critique", 33, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="dist4", 34, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="dist3", 35, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="dist2", 36, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="dist1", 37, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$names=="dist0", 38, pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$nice.names=="Ideological Distance", 39, pdata_robustness1$y)

pdata_robustness1$y <- ifelse(pdata_robustness1$specification=="clogit", pdata_robustness1$y+.25,
                              pdata_robustness1$y)
pdata_robustness1$y <- ifelse(pdata_robustness1$specification=="cjoint", pdata_robustness1$y-.25,
                              pdata_robustness1$y)

write.csv(pdata_robustness1, "pdata_robustness1.csv", row.names=FALSE) #save as CSV

#create plot
p <- ggplot(data=pdata_robustness1, group=specification)
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=.5, ymax=6.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=7.5, ymax=10.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=11.5, ymax=13.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=14.5, ymax=17.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=18.5, ymax=20.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=21.5, ymax=24.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=25.5, ymax=27.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=28.5, ymax=32.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=33.5, ymax=38.5), fill="grey")
p <- p + geom_segment(aes(x=0, xend=0, y=0, yend=38.5))
p <- p + geom_segment(aes(x=0, xend=0, y=0, yend=38.5))
p <- p + geom_vline(xintercept=-.1, color="white")
p <- p + geom_vline(xintercept=-.2, color="white")
p <- p + geom_vline(xintercept=-.3, color="white")
p <- p + geom_vline(xintercept=-.4, color="white")
p <- p + geom_vline(xintercept=-.5, color="white")
p <- p + geom_segment(aes(x=lower, xend=upper, y=y, yend=y))
p <- p + geom_point(aes(x=amce, y=y, shape=specification))
p <- p + scale_y_continuous("", position="right",
                            breaks=pdata_robustness1[specification=="clogit"&names!="header", y-.25],
                            labels=pdata_robustness1[specification=="clogit"&names!="header", nice.names])
p <- p + theme_bw()
p <- p + theme(panel.grid.major = element_blank(),
               panel.grid.minor = element_blank(),
               panel.grid.major.x = element_line(color="white"),
               panel.grid.minor.x = element_line(color="white"),
               axis.line.x = element_line(),
               panel.border = element_blank(),
               legend.position=c(.23,.085),
               legend.key.size = unit(.4, "cm"),
               text = element_text(size=14))
p <- p + scale_x_continuous("Average Marginal Component Effect",
                            breaks=seq(0, -.5, -.1))
p <- p + scale_shape_manual(name = "Regression Type",
                            breaks=c("clogit",
                                     "cjoint"),
                            labels=c("Conditional Logit",
                                     "Conjoint AMCE Estimator"),
                     values=15:16)
p <- p + annotate("text", y=pdata_robustness1[specification=="clogit"&names=="header", y-.25],
                                   x=-.25,
                  label=pdata_robustness1[specification=="clogit"&names=="header", nice.names])
p <- p + coord_cartesian(xlim = c(-.5, 0),
                         ylim = c(1.5, 37.6),
                         clip = "on")
p
ggsave("FigureA1.pdf", units="cm", width=17.7, height=17.7)

##############################################
# Figure A2: Robustness Rating Match Choices #
##############################################
#share of inconsistent choices in estimation sample (as mentioned in Online Appendix 4)
dat %>%
  filter(complete.cases(rating, chosen, dist)) %>% # restrict to our used sample for estimation including dist variable
  pivot_wider(
    id_cols = c(id_g, id_screen, screen),
    names_from = chosen,
    names_prefix = "rating_",
    values_from = rating
  ) %>%
  mutate(incorrect_order = (rating_1 < rating_0)) %>%
  count(incorrect_order) %>%
  mutate(freq = n/sum(n))
# -> same for our analysis sample

# we now re-run everything excluding respondents with any incorrect ordering (same did multiple orderings)
# extract id of incorrect orderings
id_g_incorrect_ordering <- dat %>%
  filter(complete.cases(rating, chosen)) %>%
  pivot_wider(
    id_cols = c(id_g, id_screen, screen),
    names_from = chosen,
    names_prefix = "rating_",
    values_from = rating,
    values_fn = mean
  ) %>%
  mutate(incorrect_order = (rating_1 < rating_0)) %>%
  filter(incorrect_order== TRUE) %>%
  pull(id_g)

dat_correct_ratings <- dat %>%
  filter(!id_g %in% id_g_incorrect_ordering)

# exclude missings
dat_correct_ratings.cbc <- dat_correct_ratings[complete.cases(chosen, dist), ]
length(unique(dat_correct_ratings.cbc$id_g))
# sample size (respondents) of 3395

model_clogit_correct_ratings <- clogit(chosen~
                         gender
                       +age
                       +job
                       +role
                       +critique
                       +parliament
                       +conference
                       +reform
                       +dist
                       +strata(id_screen)
                       +cluster(id_g)
                       ,data=dat_correct_ratings.cbc, method="efron", robust=TRUE)

summary(model_clogit_correct_ratings)

# custom function to get data for AMCE plot in consistent format
compute_amce_clogit_correct_ratings <- function(model) {
  coefs <- coef(model) #get coefficients
  ses <- summary(model)$coefficients[,4] #get ROBUST ses
  names <- rownames(summary(model)$coefficients) #get names of coefs
  pdata <- data.table(amce=(1/(1+exp(-coefs)))-.5, #compute probability difference to all 0's scenario
                      lower=(1/(1+exp(-(coefs-(1.96*ses)))))-.5, #compute corresponding CI
                      upper=(1/(1+exp(-(coefs+(1.96*ses)))))-.5, #compute corresponding CI
                      names=names,
                      coefs = coefs,
                      ses = ses,
                      specification = "clogit_correct_ratings"
  )
  return(pdata)
}

pdata_clogit_correct_ratings <- compute_amce_clogit_correct_ratings(model = model_clogit_correct_ratings)

pdata.model <- compute.pdata(model = model)
setnames(pdata.model, "mean", "amce")
pdata.model$specification <- "model"

pdata_robustness2 <- rbindlist(list(pdata.model, pdata_clogit_correct_ratings), fill=TRUE)

#set nice names
pdata_robustness2$nice.names <- ifelse(pdata_robustness2$names=="dist1" , "1", "")
pdata_robustness2$nice.names <- ifelse(pdata_robustness2$names=="dist2" , "2", pdata_robustness2$nice.names)
pdata_robustness2$nice.names <- ifelse(pdata_robustness2$names=="dist3" , "3", pdata_robustness2$nice.names)
pdata_robustness2$nice.names <- ifelse(pdata_robustness2$names=="dist4" , "4", pdata_robustness2$nice.names)
pdata_robustness2$nice.names <- ifelse(pdata_robustness2$names=="critiqueFormer party leader" ,
                                       "Former party leader", pdata_robustness2$nice.names)
pdata_robustness2$nice.names <- ifelse(pdata_robustness2$names=="critiqueParty factions" ,
                                       "Party factions", pdata_robustness2$nice.names)
pdata_robustness2$nice.names <- ifelse(pdata_robustness2$names=="critiqueRank-and-file members",
                                       "Rank-and-file members", pdata_robustness2$nice.names)
pdata_robustness2$nice.names <- ifelse(pdata_robustness2$names=="parliamentDivided voting",
                                       "Divided voting", pdata_robustness2$nice.names)
pdata_robustness2$nice.names <- ifelse(pdata_robustness2$names=="conferenceNeither united nor divided",
                                       "Neither united nor divided", pdata_robustness2$nice.names)
pdata_robustness2$nice.names <- ifelse(pdata_robustness2$names=="conferenceDivided" ,
                                       "Divided", pdata_robustness2$nice.names)
pdata_robustness2$nice.names <- ifelse(pdata_robustness2$names=="reformLow" , "Low", pdata_robustness2$nice.names)
pdata_robustness2$nice.names <- ifelse(pdata_robustness2$names=="roleOpposition party",
                                       "Opposition party", pdata_robustness2$nice.names)
pdata_robustness2$nice.names <- ifelse(pdata_robustness2$names=="roleJunior coalition partner",
                                       "Junior coalition partner", pdata_robustness2$nice.names)
pdata_robustness2$nice.names <- ifelse(pdata_robustness2$names=="genderMale",
                                       "Male", pdata_robustness2$nice.names)
pdata_robustness2$nice.names <- ifelse(pdata_robustness2$names=="age56 years" ,
                                       "56 years", pdata_robustness2$nice.names)
pdata_robustness2$nice.names <- ifelse(pdata_robustness2$names=="age74 years" ,
                                       "74 years", pdata_robustness2$nice.names)
pdata_robustness2$nice.names <- ifelse(pdata_robustness2$names=="jobActivist" ,
                                       "Activist", pdata_robustness2$nice.names)
pdata_robustness2$nice.names <- ifelse(pdata_robustness2$names=="jobLawyer" ,
                                       "Lawyer", pdata_robustness2$nice.names)
pdata_robustness2$nice.names <- ifelse(pdata_robustness2$names=="jobPolitician" ,
                                       "Politician", pdata_robustness2$nice.names)
pdata_robustness2$nice.names <- ifelse(pdata_robustness2$names=="jobEntrepreneur" ,
                                       "Entrepreneur", pdata_robustness2$nice.names)
pdata_robustness2$nice.names <- ifelse(pdata_robustness2$names=="jobEmployee (retired)",
                                       "Employee (retired)", pdata_robustness2$nice.names)


#addd new lines for reference categories and headers
additional.lines1 <- rbindlist(list(
                   data.table(names="dist0", nice.names="0"),
                   data.table(names="critiqueNone", nice.names="None"),
                   data.table(names="parliamentUnited", nice.names="United voting"),
                   data.table(names="conferenceUnited", nice.names="United"),
                   data.table(names="reformHigh", nice.names="High"),
                   data.table(names="rolePM party", nice.names="PM party"),
                   data.table(names="genderFemale", nice.names="Female"),
                   data.table(names="jobEmployee", nice.names="Employee"),
                   data.table(names="age38 years", nice.names="38 years"),
                   data.table(names="header", nice.names="Ideological Distance"),
                   data.table(names="header", nice.names="Party Role"),
                   data.table(names="header", nice.names="Intra-Party Critique"),
                   data.table(names="header", nice.names="Parliamentary Voting Behavior"),
                   data.table(names="header", nice.names="Behavior at Party Congress"),
                   data.table(names="header", nice.names="Reform Clarity"),
                   data.table(names="header", nice.names="Candidate: Gender"),
                   data.table(names="header", nice.names="Candidate: Age"),
                   data.table(names="header", nice.names="Candidate: Occupation")))
additional.lines1$amce <- ifelse(additional.lines1$names!="header", 0, NA)
additional.lines2 <- additional.lines1
additional.lines1$specification <- "model"
additional.lines2$specification <- "clogit_correct_ratings"
pdata_robustness2 <- rbindlist(list(pdata_robustness2, additional.lines1, additional.lines2),
                               use.names=TRUE, fill=TRUE)


## # check new names versus those from the analysis
## pdata_robustness2[,list(names, nice.names)]

#assign correct y values
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="jobActivist", 1, NA)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="jobEmployee (retired)", 2, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="jobPolitician", 3, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="jobLawyer", 4, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="jobEntrepreneur", 5, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="jobEmployee", 6, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$nice.names=="Candidate: Occupation", 7, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="age74 years", 8, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="age56 years", 9, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="age38 years", 10, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$nice.names=="Candidate: Age", 11, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="genderMale", 12, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="genderFemale", 13, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$nice.names=="Candidate: Gender", 14, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="roleOpposition party", 15, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="roleJunior coalition partner", 16, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="rolePM party", 17, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$nice.names=="Party Role", 18, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="reformLow", 19, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="reformHigh", 20, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$nice.names=="Reform Clarity", 21, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="conferenceDivided", 22, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="conferenceNeither united nor divided", 23, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="conferenceUnited", 24, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$nice.names=="Behavior at Party Congress", 25, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="parliamentDivided voting", 26, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="parliamentUnited", 27, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$nice.names=="Parliamentary Voting Behavior", 28, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="critiqueParty factions", 29, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="critiqueFormer party leader", 30, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="critiqueRank-and-file members", 31, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="critiqueNone", 32, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$nice.names=="Intra-Party Critique", 33, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="dist4", 34, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="dist3", 35, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="dist2", 36, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="dist1", 37, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$names=="dist0", 38, pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$nice.names=="Ideological Distance", 39, pdata_robustness2$y)

pdata_robustness2$y <- ifelse(pdata_robustness2$specification=="model", pdata_robustness2$y+.25,
                              pdata_robustness2$y)
pdata_robustness2$y <- ifelse(pdata_robustness2$specification=="clogit_correct_ratings", pdata_robustness2$y-.25,
                              pdata_robustness2$y)

p <- ggplot(data=pdata_robustness2, group=specification)
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=.5, ymax=6.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=7.5, ymax=10.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=11.5, ymax=13.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=14.5, ymax=17.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=18.5, ymax=20.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=21.5, ymax=24.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=25.5, ymax=27.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=28.5, ymax=32.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=33.5, ymax=38.5), fill="grey")
p <- p + geom_segment(aes(x=0, xend=0, y=0, yend=38.5))
p <- p + geom_segment(aes(x=0, xend=0, y=0, yend=38.5))
p <- p + geom_vline(xintercept=-.1, color="white")
p <- p + geom_vline(xintercept=-.2, color="white")
p <- p + geom_vline(xintercept=-.3, color="white")
p <- p + geom_vline(xintercept=-.4, color="white")
p <- p + geom_vline(xintercept=-.5, color="white")
p <- p + geom_segment(aes(x=lower, xend=upper, y=y, yend=y))
p <- p + geom_point(aes(x=amce, y=y, shape=specification))
p <- p + scale_y_continuous("", position="right",
                            breaks=pdata_robustness2[specification=="model"&names!="header", y-.25],
                            labels=pdata_robustness2[specification=="model"&names!="header", nice.names])
p <- p + theme_bw()
p <- p + theme(panel.grid.major = element_blank(),
               panel.grid.minor = element_blank(),
               panel.grid.major.x = element_line(color="white"),
               panel.grid.minor.x = element_line(color="white"),
               axis.line.x = element_line(),
               panel.border = element_blank(),
               legend.position=c(.23,.085),
               legend.key.size = unit(.4, "cm"),
               text = element_text(size=14),)
p <- p + scale_x_continuous("Average Marginal Component Effect",
                            breaks=seq(0, -.5, -.1))
p <- p + scale_shape_manual(name = "Sample type",
                            breaks=c("model",
                                     "clogit_correct_ratings"),
                            labels=c("Full sample",
                                     "Consistent ratings"),
                     values=15:16)
p <- p + annotate("text", y=pdata_robustness2[specification=="model"&names=="header", y-.25],
                                   x=-.25,
                  label=pdata_robustness2[specification=="model"&names=="header", nice.names])
p <- p + coord_cartesian(xlim = c(-.5, 0),
                         ylim = c(1.5, 37.7),
                         clip = "on")
p
ggsave("FigureA2.pdf", units="cm", width=17.7, height=17.7)


###################################################################
# Figure A3: Only with choices who saw all attributes of critique #
###################################################################

 # find attributes of critique respondent has seen
 dat <- dat[,c1:=ifelse(screen>=2, critique[party==1], NA), id_g]
 dat <- dat[,c2:=ifelse(screen>=2, critique[party==2], NA), id_g]
 dat <- dat[,c3:=ifelse(screen>=3, critique[party==3], NA), id_g]
 dat <- dat[,c4:=ifelse(screen>=3, critique[party==4], NA), id_g]
 dat <- dat[,c5:=ifelse(screen>=4, critique[party==5], NA), id_g]
 dat <- dat[,c6:=ifelse(screen>=4, critique[party==6], NA), id_g]
 dat <- dat[,c7:=ifelse(screen>=5, critique[party==7], NA), id_g]
 dat <- dat[,c8:=ifelse(screen>=5, critique[party==8], NA), id_g]
 dat <- dat[,c9:=ifelse(screen>=6, critique[party==9], NA), id_g]
 dat <- dat[,c10:=ifelse(screen>=6, critique[party==10], NA), id_g]
 dat <- dat[,c11:=ifelse(screen>=7, critique[party==11], NA), id_g]
 dat <- dat[,c12:=ifelse(screen>=7, critique[party==12], NA), id_g]
 dat <- dat[,c13:=ifelse(screen>=8, critique[party==12], NA), id_g]
 dat <- dat[,c14:=ifelse(screen>=8, critique[party==14], NA), id_g]
 dat <- dat[,c15:=ifelse(screen>=9, critique[party==15], NA), id_g]
 dat <- dat[,c16:=ifelse(screen>=9, critique[party==16], NA), id_g]
 dat <- dat[,c17:=ifelse(screen>=10, critique[party==17], NA), id_g]
 dat <- dat[,c18:=ifelse(screen>=10, critique[party==18], NA), id_g]

#unique number of critique attributes repsondent has seen
dat <- dat[,critique.seen:=length(unique(na.omit(c(c1, c2, c3, c4, c5, c6, c7, c8, c9, c10,
                      c11, c12, c13, c14, c15, c16, c17, c18)))), .(id_g, party)]

#remaining choices
table(dat$critique.seen==4)[2]/2

table(dat$critique.seen==4)[2]/sum(table(dat$critique.seen==4))

model.seen <- clogit(chosen~
                  dist +
                  critique +
                  parliament +
                  conference +
                  reform +
                  role +
                  gender +
                  age +
                  job +
                  strata(id_screen) +
                  cluster(id_g),
                data=dat[critique.seen==4], method="efron", robust=TRUE)
summary(model.seen)

compute.pdata <- function(model) {
    #(get data for) AMCEs plot
    coefs <- coef(model) #get cofficients
    ses <- summary(model)$coefficients[,4] #get ROBUST res
    names <- rownames(summary(model)$coefficients) #get names of coefs
    pdata <- data.table(mean=(1/(1+exp(-coefs)))-.5, #compute probability difference to all 0's scenario
                        lower=(1/(1+exp(-(coefs-(1.96*ses)))))-.5, #compute corresponding CI
                        upper=(1/(1+exp(-(coefs+(1.96*ses)))))-.5, #compute corresponding CI
                        names=names
                        )
    return(pdata)
}

pdata.seen <- compute.pdata(model = model.seen)
pdata.seen$model <- "model.seen"
pdata.model <- compute.pdata(model = model)
pdata.model$model <- "model"

pdata <- rbindlist(list(pdata.seen, pdata.model))

#set nice names
pdata$nice.names <- ifelse(pdata$names=="dist1" , "1", "")
pdata$nice.names <- ifelse(pdata$names=="dist2" , "2", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="dist3" , "3", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="dist4" , "4", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="critiqueFormer party leader" ,
                                       "Former party leader", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="critiqueParty factions" ,
                                       "Party factions", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="critiqueRank-and-file members",
                                       "Rank-and-file members", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="parliamentDivided voting",
                                       "Divided voting", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="conferenceNeither united nor divided",
                                       "Neither united nor divided", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="conferenceDivided" ,
                                       "Divided", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="reformLow" , "Low", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="roleOpposition party",
                                       "Opposition party", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="roleJunior coalition partner",
                                       "Junior coalition partner", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="genderMale",
                                       "Male", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="age56 years" ,
                                       "56 years", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="age74 years" ,
                                       "74 years", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="jobActivist" ,
                                       "Activist", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="jobLawyer" ,
                                       "Lawyer", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="jobPolitician" ,
                                       "Politician", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="jobEntrepreneur" ,
                                       "Entrepreneur", pdata$nice.names)
pdata$nice.names <- ifelse(pdata$names=="jobEmployee (retired)",
                                       "Employee (retired)", pdata$nice.names)

#addd new lines for reference categories and headers
additional.lines1 <- rbindlist(list(
                   data.table(mean=0, names="dist0", nice.names="0"),
                   data.table(mean=0, names="critiqueNone", nice.names="None"),
                   data.table(mean=0, names="parliamentUnited", nice.names="United voting"),
                   data.table(mean=0, names="conferenceUnited", nice.names="United"),
                   data.table(mean=0, names="reformHigh", nice.names="High"),
                   data.table(mean=0, names="rolePM party", nice.names="PM party"),
                   data.table(mean=0, names="genderFemale", nice.names="Female"),
                   data.table(mean=0, names="jobEmployee", nice.names="Employee"),
                   data.table(mean=0, names="age38 years", nice.names="38 years"),
                   data.table(mean=NA, names="header", nice.names="Ideological Distance"),
                   data.table(mean=NA, names="header", nice.names="Party Role"),
                   data.table(mean=NA, names="header", nice.names="Intra-Party Critique"),
                   data.table(mean=NA, names="header", nice.names="Parliamentary Voting Behavior"),
                   data.table(mean=NA, names="header", nice.names="Behavior at Party Congress"),
                   data.table(mean=NA, names="header", nice.names="Reform Clarity"),
                   data.table(mean=NA, names="header", nice.names="Candidate: Gender"),
                   data.table(mean=NA, names="header", nice.names="Candidate: Age"),
                   data.table(mean=NA, names="header", nice.names="Candidate: Occupation")))
additional.lines2 <- additional.lines1
additional.lines1$model <- "model"
additional.lines2$model <- "model.seen"
pdata <- rbindlist(list(pdata, additional.lines1, additional.lines2),
                               use.names=TRUE, fill=TRUE)


## # check new names versus those from the analysis
## pdata[,list(names, nice.names)]

#assign correct y values
pdata$y <- ifelse(pdata$names=="jobActivist", 1, NA)
pdata$y <- ifelse(pdata$names=="jobEmployee (retired)", 2, pdata$y)
pdata$y <- ifelse(pdata$names=="jobPolitician", 3, pdata$y)
pdata$y <- ifelse(pdata$names=="jobLawyer", 4, pdata$y)
pdata$y <- ifelse(pdata$names=="jobEntrepreneur", 5, pdata$y)
pdata$y <- ifelse(pdata$names=="jobEmployee", 6, pdata$y)
pdata$y <- ifelse(pdata$nice.names=="Candidate: Occupation", 7, pdata$y)
pdata$y <- ifelse(pdata$names=="age74 years", 8, pdata$y)
pdata$y <- ifelse(pdata$names=="age56 years", 9, pdata$y)
pdata$y <- ifelse(pdata$names=="age38 years", 10, pdata$y)
pdata$y <- ifelse(pdata$nice.names=="Candidate: Age", 11, pdata$y)
pdata$y <- ifelse(pdata$names=="genderMale", 12, pdata$y)
pdata$y <- ifelse(pdata$names=="genderFemale", 13, pdata$y)
pdata$y <- ifelse(pdata$nice.names=="Candidate: Gender", 14, pdata$y)
pdata$y <- ifelse(pdata$names=="roleOpposition party", 15, pdata$y)
pdata$y <- ifelse(pdata$names=="roleJunior coalition partner", 16, pdata$y)
pdata$y <- ifelse(pdata$names=="rolePM party", 17, pdata$y)
pdata$y <- ifelse(pdata$nice.names=="Party Role", 18, pdata$y)
pdata$y <- ifelse(pdata$names=="reformLow", 19, pdata$y)
pdata$y <- ifelse(pdata$names=="reformHigh", 20, pdata$y)
pdata$y <- ifelse(pdata$nice.names=="Reform Clarity", 21, pdata$y)
pdata$y <- ifelse(pdata$names=="conferenceDivided", 22, pdata$y)
pdata$y <- ifelse(pdata$names=="conferenceNeither united nor divided", 23, pdata$y)
pdata$y <- ifelse(pdata$names=="conferenceUnited", 24, pdata$y)
pdata$y <- ifelse(pdata$nice.names=="Behavior at Party Congress", 25, pdata$y)
pdata$y <- ifelse(pdata$names=="parliamentDivided voting", 26, pdata$y)
pdata$y <- ifelse(pdata$names=="parliamentUnited", 27, pdata$y)
pdata$y <- ifelse(pdata$nice.names=="Parliamentary Voting Behavior", 28, pdata$y)
pdata$y <- ifelse(pdata$names=="critiqueParty factions", 29, pdata$y)
pdata$y <- ifelse(pdata$names=="critiqueFormer party leader", 30, pdata$y)
pdata$y <- ifelse(pdata$names=="critiqueRank-and-file members", 31, pdata$y)
pdata$y <- ifelse(pdata$names=="critiqueNone", 32, pdata$y)
pdata$y <- ifelse(pdata$nice.names=="Intra-Party Critique", 33, pdata$y)
pdata$y <- ifelse(pdata$names=="dist4", 34, pdata$y)
pdata$y <- ifelse(pdata$names=="dist3", 35, pdata$y)
pdata$y <- ifelse(pdata$names=="dist2", 36, pdata$y)
pdata$y <- ifelse(pdata$names=="dist1", 37, pdata$y)
pdata$y <- ifelse(pdata$names=="dist0", 38, pdata$y)
pdata$y <- ifelse(pdata$nice.names=="Ideological Distance", 39, pdata$y)

pdata$y <- ifelse(pdata$model=="model", pdata$y+.25,
                              pdata$y)
pdata$y <- ifelse(pdata$model=="model.seen", pdata$y-.25,
                              pdata$y)

p <- ggplot(data=pdata, group=model)
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=.5, ymax=6.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=7.5, ymax=10.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=11.5, ymax=13.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=14.5, ymax=17.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=18.5, ymax=20.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=21.5, ymax=24.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=25.5, ymax=27.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=28.5, ymax=32.5), fill="grey")
p <- p + geom_rect(aes(xmin=.02, xmax=-.55, ymin=33.5, ymax=38.5), fill="grey")
p <- p + geom_segment(aes(x=0, xend=0, y=0, yend=38.5))
p <- p + geom_segment(aes(x=0, xend=0, y=0, yend=38.5))
p <- p + geom_vline(xintercept=-.1, color="white")
p <- p + geom_vline(xintercept=-.2, color="white")
p <- p + geom_vline(xintercept=-.3, color="white")
p <- p + geom_vline(xintercept=-.4, color="white")
p <- p + geom_vline(xintercept=-.5, color="white")
p <- p + geom_segment(aes(x=lower, xend=upper, y=y, yend=y))
p <- p + geom_point(aes(x=mean, y=y, shape=model))
p <- p + scale_y_continuous("", position="right",
                            breaks=pdata[model=="model"&names!="header", y-.25],
                            labels=pdata[model=="model"&names!="header", nice.names])
p <- p + theme_bw()
p <- p + theme(panel.grid.major = element_blank(),
               panel.grid.minor = element_blank(),
               panel.grid.major.x = element_line(color="white"),
               panel.grid.minor.x = element_line(color="white"),
               axis.line.x = element_line(),
               panel.border = element_blank(),
               legend.position=c(.2,.095),
               legend.key.size = unit(.4, "cm"),
               text = element_text(size=14))
p <- p + scale_x_continuous("Average Marginal Component Effect",
                            breaks=seq(0, -.5, -.1))
p <- p + scale_shape_manual(name = "Intra-party critique\nattributes seen",
                            breaks=c("model",
                                     "model.seen"),
                            labels=c("Any numer",
                                     "All"),
                     values=15:16)
p <- p + annotate("text", y=pdata[model=="model"&names=="header", y-.25],
                                   x=-.25,
                  label=pdata[model=="model"&names=="header", nice.names])
p <- p + coord_cartesian(xlim = c(-.48, -.01),
                         ylim = c(1.5, 37.6),
                         clip = "on")
p
ggsave("FigureA3.pdf", units="cm", width=17.7, height=17.7)

#############
# Figure A4 #
#############
#compute predicted probabilities for all subplots
#set scenarios
scenarios=expand.grid(role=c("PM party"),
                          conference= c("United", "Neither united nor divided", "Divided"),
                          parliament=c("United voting", "Divided voting"),
                          critique= c("Rank-and-file members", "Former party leader", "Party factions", "None"),
                          reform=c("High"),
                          gender=c("Female"),
                          age=c("38 years"),
                          job=c( "Employee"),
                          dist=c("0", "1", "2", "3", "4"),
                          id_screen=1)
scenarios <- as.data.table(scenarios)

# create design matrix that has indicator vars instead of factors
scenarios.expanded <- model.matrix(model, data=scenarios)

# compute probabilities
linpred <- scenarios.expanded%*%t(betas)
probs <- exp(linpred)/(1+exp(linpred))

# find median probability and CIs for plotting
scenarios$prob <- apply(probs, 1, median)*100
scenarios$upper <- apply(probs, 1, quantile, p=.975)*100
scenarios$lower <- apply(probs, 1, quantile, p=.025)*100


# a) Effect of ideological distance by unity in parliamentary voting
pdata <- scenarios[role == "PM party" &
                   conference == "United" &
                   parliament %in% c("United voting", "Divided voting") &
                   critique == "None" &
                   reform == "High" &
                   gender == "Female" &
                   age == "38 years" &
                   job == "Employee" &
                   dist %in% c(0:4),
                   .(parliament, dist, lower, prob, upper)]
pdata$dist <- as.numeric(as.character(pdata$dist))
pdata$dist <- ifelse(pdata$parliament=="United voting", pdata$dist-.1, pdata$dist)
pdata$dist <- ifelse(pdata$parliament=="Divided voting", pdata$dist+.1, pdata$dist)

p <- ggplot(data=pdata, aes(x=dist))
p <- p + geom_rect(aes(xmin=-.6, xmax=.4, ymin=0, ymax=55), fill="grey")
p <- p + geom_rect(aes(xmin=.6, xmax=1.4, ymin=0, ymax=55), fill="grey")
p <- p + geom_rect(aes(xmin=1.6, xmax=2.4, ymin=0, ymax=55), fill="grey")
p <- p + geom_rect(aes(xmin=2.6, xmax=3.4, ymin=0, ymax=55), fill="grey")
p <- p + geom_rect(aes(xmin=3.6, xmax=4.4, ymin=0, ymax=55), fill="grey")
p <- p + geom_hline(yintercept=10, color="white")
p <- p + geom_hline(yintercept=20, color="white")
p <- p + geom_hline(yintercept=30, color="white")
p <- p + geom_hline(yintercept=40, color="white")
p <- p + geom_hline(yintercept=50, color="white")
p <- p + geom_segment(aes(x=dist, xend=dist, y=lower, yend=upper), size=1)
p <- p + geom_point(aes(y=prob, shape=parliament), size=3)
p <- p + xlab("Ideological Distance")
p <- p + scale_y_continuous("Predicted Vote Share\nvs. Baseline",
                            breaks=seq(0, 50, 10))
p <- p + scale_shape_discrete(name = "Voting Behavior in Parliament",
                              labels=c("United", "Divided"))
p <- p + theme_bw()
p <- p + theme(panel.grid.major.x = element_blank(),
               panel.grid.minor.x = element_blank(),
               axis.line = element_line(),
               panel.border = element_blank(),
               legend.position=c(.859,.89),
               legend.background = element_rect(linetype="solid", color="black"),
               text = element_text(size=14))
p <- p + coord_cartesian(xlim = c(-.5, 4.3),
                         ylim = c(1, 50),
                         clip = "on")
p
ggsave("FigureA4_1.pdf", width=32.3, height=12.1, units="cm")


# b) Effect of ideological distance by unity in party conference behavior
pdata <- scenarios[role == "PM party" &
                         conference %in% c("United", "Neither united nor divided", "Divided") &
                         parliament == "United voting" &
                         critique == "None" &
                         reform == "High" &
                         gender == "Female" &
                         age == "38 years" &
                         job == "Employee" &
                         dist %in% c(0:4),
                         .(conference, dist, lower, prob, upper)]
pdata$dist <- as.numeric(as.character(pdata$dist))
pdata$dist <- ifelse(pdata$conference=="United", pdata$dist-.1, pdata$dist)
pdata$dist <- ifelse(pdata$conference=="Divided", pdata$dist+.1, pdata$dist)

p <- ggplot(data=pdata, aes(x=dist))
p <- p + geom_rect(aes(xmin=-.6, xmax=.4, ymin=0, ymax=55), fill="grey")
p <- p + geom_rect(aes(xmin=.6, xmax=1.4, ymin=0, ymax=55), fill="grey")
p <- p + geom_rect(aes(xmin=1.6, xmax=2.4, ymin=0, ymax=55), fill="grey")
p <- p + geom_rect(aes(xmin=2.6, xmax=3.4, ymin=0, ymax=55), fill="grey")
p <- p + geom_rect(aes(xmin=3.6, xmax=4.4, ymin=0, ymax=55), fill="grey")
p <- p + geom_hline(yintercept=10, color="white")
p <- p + geom_hline(yintercept=20, color="white")
p <- p + geom_hline(yintercept=30, color="white")
p <- p + geom_hline(yintercept=40, color="white")
p <- p + geom_hline(yintercept=50, color="white")
p <- p + geom_segment(aes(x=dist, xend=dist, y=lower, yend=upper), size=1)
p <- p + geom_point(aes(y=prob, shape=conference), size=3)
p <- p + xlab("Ideological Distance")
p <- p + scale_y_continuous("Predicted Vote Share\nvs. Baseline",
                            breaks=seq(0, 50, 10))
p <- p + scale_shape_discrete(name = "Behavior at Party Congress")
p <- p + theme_bw()
p <- p + theme(panel.grid.major.x = element_blank(),
               panel.grid.minor.x = element_blank(),
               axis.line = element_line(),
               panel.border = element_blank(),
               legend.position=c(.866,.86),
               legend.background = element_rect(linetype="solid", color="black"),
               text = element_text(size=14))
p <- p + coord_cartesian(xlim = c(-.5, 4.3),
                         ylim = c(1, 50),
                         clip = "on")
p
ggsave("FigureA4_2.pdf", width=32.3, height=12.1, units="cm")

# c) Effect of ideological distance by internal critique
pdata <- scenarios[role == "PM party" &
                                conference == "United" &
                              parliament == "United voting" &
                              critique %in% c("Rank-and-file members"
                                              ,"Former party leader"
                                              ,"Party factions"
                                              ,"None") &
                         reform == "High" &
                         gender == "Female" &
                         age == "38 years" &
                         job == "Employee" &
                         dist %in% c(0:4),
                         .(critique, dist, lower, prob, upper)]
pdata$dist <- as.numeric(as.character(pdata$dist))
pdata$dist <- ifelse(pdata$critique=="Rank-and-file members", pdata$dist-.2, pdata$dist)
pdata$dist <- ifelse(pdata$critique=="None", pdata$dist-.1, pdata$dist)
pdata$dist <- ifelse(pdata$critique=="Former leader", pdata$dist+.1, pdata$dist)
pdata$dist <- ifelse(pdata$critique=="Party factions", pdata$dist+.2, pdata$dist)


p <- ggplot(data=pdata, aes(x=dist))
p <- p + geom_rect(aes(xmin=-.6, xmax=.4, ymin=0, ymax=55), fill="grey")
p <- p + geom_rect(aes(xmin=.6, xmax=1.4, ymin=0, ymax=55), fill="grey")
p <- p + geom_rect(aes(xmin=1.6, xmax=2.4, ymin=0, ymax=55), fill="grey")
p <- p + geom_rect(aes(xmin=2.6, xmax=3.4, ymin=0, ymax=55), fill="grey")
p <- p + geom_rect(aes(xmin=3.6, xmax=4.4, ymin=0, ymax=55), fill="grey")
p <- p + geom_hline(yintercept=10, color="white")
p <- p + geom_hline(yintercept=20, color="white")
p <- p + geom_hline(yintercept=30, color="white")
p <- p + geom_hline(yintercept=40, color="white")
p <- p + geom_hline(yintercept=50, color="white")
p <- p + geom_segment(aes(x=dist, xend=dist, y=lower, yend=upper), size=1)
p <- p + geom_point(aes(y=prob, shape=critique), size=3)
p <- p + xlab("Ideological Distance")
p <- p + scale_y_continuous("Predicted Vote Share\nvs. Baseline",
                            breaks=seq(0, 50, 10))
p <- p + scale_shape_manual(name = "Intra-Party Critique",
                            breaks=c("Rank-and-file members",
                                     "None",
                                     "Former party leader",
                                     "Party factions"),
                     values=15:18)
p <- p + theme_bw()
p <- p + theme(panel.grid.major.x = element_blank(),
               panel.grid.minor.x = element_blank(),
               axis.line = element_line(),
               panel.border = element_blank(),
               legend.position=c(.886,.832),
               legend.background = element_rect(linetype="solid", color="black"),
               text = element_text(size=14))
p <- p + coord_cartesian(xlim = c(-.5, 4.3),
                         ylim = c(1, 50),
                         clip = "on")
p
ggsave("FigureA4_3.pdf", width=32.3, height=12.1, units="cm")

#########################################
# Information on Software versions used #
#########################################
print("Information on Software versions used:")

sessionInfo()

closeAllConnections() #close log file

